program brentry_bioprobit_neg
version 10.1
args lnf s1 v1 f1 alpha1_1 gamma1_1 s2 v2 f2 alpha1_2 alpha2_2 alpha3_2 gamma1_2 gamma2_2 gamma3_2 r

tempvar p1u p2u p3u p1h p11 p12 p13 pi_h1 pi_u1 pi_u2 pi_u3
tempname rho

qui gen double `rho' = tanh(`r')

qui gen double `pi_h1'=`s1'*(`v1'+`alpha1_1')-`f1'-`gamma1_1'

qui gen double `pi_u1'=`s2'*(`v2'+`alpha1_2')-`f2'-`gamma1_2'
qui gen double `pi_u2'=`s2'*(`v2'+`alpha1_2'-`alpha2_2')-`f2'-`gamma1_2'-`gamma2_2'
qui gen double `pi_u3'=`s2'*(`v2'+`alpha1_2'-`alpha2_2'-`alpha3_2')-`f2'-`gamma1_2'-`gamma2_2'-`gamma3_2'

qui gen double `p1u'=normal(-1*`pi_u1')
qui gen double `p2u'=normal(-1*`pi_u2')
qui gen double `p3u'=normal(-1*`pi_u3')

qui gen double `p1h'=normal(-1*`pi_h1')	
qui gen double `p11'=binormal(-1*`pi_h1', -1*`pi_u1', `rho')
qui gen double `p12'=binormal(-1*`pi_h1', -1*`pi_u2', `rho')
qui gen double `p13'=binormal(-1*`pi_h1', -1*`pi_u3', `rho')

qui replace `lnf' = ln(`p11') if $ML_y1==0 & $ML_y4==0
qui replace `lnf' = ln(`p12'-`p11') if $ML_y1==0 & $ML_y4==1
qui replace `lnf' = ln(`p13'-`p12') if $ML_y1==0 & $ML_y4==2
qui replace `lnf' = ln(`p1h'-`p13') if $ML_y1==0 & $ML_y4>=3

qui replace `lnf' = ln(`p1u'-`p11') if $ML_y1>=1 & $ML_y4==0
qui replace `lnf' = ln(`p2u'-`p1u'-`p12'+`p11') if $ML_y1>=1 & $ML_y4==1
qui replace `lnf' = ln(`p3u'-`p2u'-`p13'+`p12') if $ML_y1>=1 & $ML_y4==2
qui replace `lnf' = ln(1-`p1h'-`p3u'+`p13') if $ML_y1>=1 & $ML_y4>=3


end

